/*=======================================================================================
	Dynamics.do
	
		Uses the Personal Finance Website (PFW) data to analyze the response of 
		consumption to the APFD payments
		
	Author: Lorenz Kueng, Sep 2017
=========================================================================================*/

cap log close PFW_04b
log using "$homedir/log-files/PFW_04b_$date.log", text replace name(PFW_04b)
	


*=============================================================
* Parametric Analysis: Figures
*=============================================================

use "$homedir/data/stata/PFW_monthly.dta", clear

	generate date = Month

	
*** Create additional variables
	
	** Controls
	
	foreach var in residentialstatus marriagestatus {
		display _n(1) "`var'"
		cap drop `var'_
		encode   `var', generate(`var'_)
		bysort userid (date): egen _temp = mode(`var'_), maxmode
		replace  `var'_ = _temp if `var'_==. & _temp!=. 
		replace  `var'_ = 999   if `var'_==. 
		drop _temp
	}
	foreach var in age edu profession homezipcode {
		display _n(1) "`var'"
		cap drop `var'_
		generate `var'_ = `var'
		bysort userid (date): egen _temp = mode(`var'_), maxmode
		replace  `var'_ = _temp if `var'_==. & _temp!=. 
		replace  `var'_ = 999   if `var'_==. 
		drop _temp
	}
	
	foreach var in incomeM incomeY liquid_asset1 {
		winsor `var', g(`var'W) p(0.01)
		sum `var' `var'W, d
	}

	xtset userid date

	gen D_incomeMW = D.incomeMW
	


foreach depvar of varlist nondur_serv durable { 
 
	global depvar = "`depvar'"
		
	xtset userid date
		

	*-----------------------------------------------------------
	* Winsorize dependent variable at 1%
	*-----------------------------------------------------------
			
		* save a copy of the original variable

		cap drop ${depvar}_original
		generate ${depvar}_original = ${depvar}
		
		* winsorize at 1% (conditional on non-zero observations for disaggregated categories)
		
		summarize ${depvar}, detail
		summarize ${depvar} if ${depvar}!=0, detail

		count if ${depvar}==0 
		
		if `r(N)'<5000 {
		
			cap drop temp
			winsor ${depvar}, generate(temp) p(0.01) 
			replace temp = 0 if temp==.
			replace ${depvar} = temp
		}
		else { // winsorize conditional on non-zeros for disaggregated variables
		
			cap drop temp
			winsor ${depvar} if ${depvar}!=0, generate(temp) p(0.01) 
			replace temp = 0 if temp==.
			replace ${depvar} = temp
		}
		
		* make all variables positive (i.e., expenditures)
		
		replace $depvar = abs(${depvar}) 
		summarize ${depvar}, detail


	*-----------------------------------------------------------
	* Figures: Dyanmic effects
	*-----------------------------------------------------------
	
	foreach spec in with_controls1 no_controls with_controls2 with_controls3 {
	
		lab var APFD ""

		*** Without additional controls
		
		if "`spec'"=="no_controls" { // only main effects: family size, time and state FEs
		
			reghdfe D.${depvar} L(-6/8).APFD ///
				if APFDid_annual!=. ///
				, absorb(date Alaska family_sizeImp) vce(cluster userid)
		}
		
		
		*** Adding more controls 
		
		if "`spec'"=="with_controls1" { // liquid wealth and income changes 
					
			reghdfe D.${depvar} L(-6/8).APFD ///
					D_incomeMW ///
					liquid_asset1W ///
				if APFDid_annual!=. /// NOTE: APFDid_annual!=. conditions on Alaskans for whom we know they received a dividend via direct deposit in early October
				, absorb(date Alaska family_sizeImp) vce(cluster userid) 
		}


		if "`spec'"=="with_controls2" { // liquid wealth and leads & lags of income changes 
		
			reghdfe D.${depvar} L(-6/8).APFD ///
					L(-6/8).D_incomeMW /// 
					liquid_asset1W ///
				if APFDid_annual!=. ///
				, absorb(date Alaska family_sizeImp) vce(cluster userid)
		}
		

		if "`spec'"=="with_controls3" { // liquid wealth and leads & lags of income changes and state-time FE
		
			reghdfe D.${depvar} L(-6/8).APFD ///
					L(-6/8).D_incomeMW /// 
					liquid_asset1W ///
					i.date##i.Alaska ///
				if APFDid_annual!=. ///
				, absorb(date Alaska family_sizeImp) vce(cluster userid)				
		}

		
		*-------------------
		* Plot regression coefficients
		*-------------------

		if "`spec'"=="with_controls1" { // liquid wealth and income changes 

		coefplot, keep(*APFD*) vertical ///
					coefl(	///
							F6.APFD = "-6" ///
							F5.APFD = "-5" ///
							F4.APFD = "-4" ///
							F3.APFD = "-3" ///
							F2.APFD = "-2" ///
							 F.APFD = "-1" ///
							   APFD = "0" ///
							 L.APFD = "1" ///
							L2.APFD = "2" ///
							L3.APFD = "3" ///
							L4.APFD = "4" ///
							L5.APFD = "5" ///
							L6.APFD = "6" ///
							L7.APFD = "7" ///
							L8.APFD = "8" ///
					) ///
					xtitle("months since dividend payment (event time)") ytitle("marginal effects ({&beta}s)") ///
					ciopts(recast(rcap) lc(black)) mcolor(black) mlabel format(%-4.2f) mlabc(black) ///
					yline(0, lp(solid) lc(black) lw(thin))  ///
					xline(7, lp(dash) lc(blue) lw(vthin)) /// 
					graphregion(color(white) lstyle(none)) legend(off) 

			local file "$homedir/results/figures/PFW_D${depvar}_dynamics_monthly_`spec'"
				graph export "`file'.eps", replace
				cap rm       "`file'.pdf"
				!epstopdf    "`file'.eps"
				rm           "`file'.eps"
				graph export "`file'.tif", replace
		
		lincom F.APFD, level(90) // for 95%-confidence upper bound on anticipation effect

		}
		
		
		*-------------------
		* Cumulative MPC
		*-------------------

		global LLL = 6 // horizon for dynamic plot (LLL=5 equals 2 quarters. Show 6, which is first months where cumulative effect is not significant anymore)
		
		capture noisily {
		
		foreach var in net_b net_se net_p cumu_b cumu_se cumu_p horizon {
			cap drop `var'
			generate `var'=.
		}
		
		local LLL = ${LLL}
		forvalues l=0(1)`LLL' {
			
			local L0 = `l'+1
			local L1 = `l'
			local L2 = `l'-1
			local L3 = `l'-2
			local L4 = `l'-3
			local L5 = `l'-4
			local L6 = `l'-5
			local L7 = `l'-6
			local L8 = `l'-7	
			forvalues i = 0(1)${LLL} {
				
				if `L`i''<0 {
					local L`i' = 0
				} // end "if"
			} // enf "forvalues i"
			
			local obs = `l'+1
			replace horizon = `l' in `obs'
			
			if `l'==0 & `l'<=${LLL} {
			capture noisily{
			
				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] in `obs'
				test 	          `L0'*   APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'
				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD] in `obs'
				test 	                  APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'	
			}
			}
			if `l'==1 & `l'<=${LLL} {
			capture noisily{
			
				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==2 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==3 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==4 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] + `L4'*_b[L4.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  + `L4'*   L4.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] +      _b[L4.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  +         L4.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==5 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] + `L4'*_b[L4.APFD] + `L5'*_b[L5.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  + `L4'*   L4.APFD  + `L5'*   L5.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] +      _b[L4.APFD] +      _b[L5.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  +         L4.APFD  +         L5.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==6 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] + `L4'*_b[L4.APFD] + `L5'*_b[L5.APFD] + `L6'*_b[L6.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  + `L4'*   L4.APFD  + `L5'*   L5.APFD  + `L6'*   L6.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] +      _b[L4.APFD] +      _b[L5.APFD] +      _b[L6.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  +         L4.APFD  +         L5.APFD  +         L6.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==7 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] + `L4'*_b[L4.APFD] + `L5'*_b[L5.APFD] + `L6'*_b[L6.APFD]  + `L7'*_b[L7.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  + `L4'*   L4.APFD  + `L5'*   L5.APFD  + `L6'*   L6.APFD   + `L7'*   L7.APFD  = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] +      _b[L4.APFD] +      _b[L5.APFD] +      _b[L6.APFD]  +      _b[L7.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  +         L4.APFD  +         L5.APFD  +         L6.APFD   +         L7.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			}
			}
			if `l'==8 & `l'<=${LLL} {
			capture noisily{

				* cumulative effect
				
				replace cumu_b  = `L0'*_b[APFD] + `L1'*_b[L1.APFD] + `L2'*_b[L2.APFD] + `L3'*_b[L3.APFD] + `L4'*_b[L4.APFD] + `L5'*_b[L5.APFD] + `L6'*_b[L6.APFD]  + `L7'*_b[L7.APFD] + `L8'*_b[L8.APFD] in `obs'
				test 	          `L0'*   APFD  + `L1'*   L1.APFD  + `L2'*   L2.APFD  + `L3'*   L3.APFD  + `L4'*   L4.APFD  + `L5'*   L5.APFD  + `L6'*   L6.APFD   + `L7'*   L7.APFD  + `L8'*   L8.APFD = 0
				replace cumu_se = cumu_b[`obs']/sqrt(r(F)) in `obs'
				replace cumu_p  = r(p)                     in `obs'

				
				* net effect: C_{t+l}-C_{t-1}
				
				replace net_b   =      _b[APFD]  +      _b[L1.APFD] +      _b[L2.APFD] +      _b[L3.APFD] +      _b[L4.APFD] +      _b[L5.APFD] +      _b[L6.APFD]  +      _b[L7.APFD] +      _b[L8.APFD] in `obs'
				test 	                  APFD   +         L1.APFD  +         L2.APFD  +         L3.APFD  +         L4.APFD  +         L5.APFD  +         L6.APFD   +         L7.APFD  +         L8.APFD  = 0 
				replace net_se  = net_b[`obs']/sqrt(r(F))  in `obs'
				replace net_p   = r(p)                     in `obs'		
			} // end "if `l'==8"
			}
		} // end "forvalues l=0(1)${LLL} "
		replace cumu_se = abs(cumu_se)
		replace net_se  = abs(net_se)
		}
				
		** Cummulative MPC
		
		cap drop cumu_2seup
		generate cumu_2seup  = cumu_b + 2*cumu_se
		cap drop cumu_2selow
		generate cumu_2selow = cumu_b - 2*cumu_se
		cap drop mlab 
		generate mlab = string(cumu_b,"%3.2f")

		twoway connected cumu_b      horizon, lp(solid) lc(black) connect(line) mc(black) mlabel(mlab) mlabc(black) mlabposition(12) || ///
			   line      cumu_2seup  horizon, lp(dash)  lc(black) || ///
			   line      cumu_2selow horizon, lp(dash)  lc(black) || ///
			, xtitle("horizon  (months)") xlabel(0(1)${LLL}) ytitle("cumulative effect") graphregion(color(white)) legend( off ) yline(0, lc(black) lw(thin)  )

			local file "$homedir/results/figures/PFW_D${depvar}_cumu_monthly_`spec'"
				graph export "`file'.eps", replace
				cap rm       "`file'.pdf"
				!epstopdf    "`file'.eps"
				rm           "`file'.eps"
				graph export "`file'.tif", replace

	} // end of spec	
} // end of depvar	

log close PFW_04b
